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The spin-dependent corrections to the static interquark potential are relevant to describing the 
fine and hyper-fine splittings of the heavy quarkonium spectra. We investigate these corrections 
in SU(3) lattice gauge theory with the Polyakov loop correlation function as the quark source by 
applying the multi-level algorithm. We observe remarkably clean signals for the spin-dependent 
potentials up to intermediate distances. 
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1. Introduction 

The spin-dependent (spin-orbit, spin-spin) interquark potentials are relevant to describing the 
fine and hyper-fine splittings of the heavy quarkonium spectra and thus it is interesting to determine 
their behavior directly from QCD. Eichten and Feinberg ^ derived in this context the general form 
of the potential including the spin-dependent corrections up to 0(1 /m 2 ), 

v{r) _ m M^\(m +2 m\J^ + i^\m 

\ 2mf 2/«2 / \ t r J \lm\m2 2m\m.2 J r 

1 f {si-r)(s 2 -r) Si-s 2 \ Si-s 2 

+ 2 r~ Mn , (i-i) 

m\ni2 V ' -J / 5m\ni2 

where the locations of quark and antiquark are r\, f 2 (r = \r\ — f 2 \). m\, ni2 (= m) denote their 
masses, s\, s 2 the spins, l\ = —l 2 = I the orbital angular momenta. Vo(r) is the spin-independent 
static potential and V, (r) (i = 1, ... ,4) the spin-dependent potentials, which are expressed in terms 
of the correlation function of two field strength operators attached to the quark and antiquark, 
respectively. The prime denotes the derivative with respect to r. 

The determination of these potentials through lattice Monte Carlo simulations goes back to the 
1980s [|2[ ||, ||, |p. The latest investigations are found in refs. f7|. The qualitative (quantitative 
to some extent) findings which seem to be established are that while the spin-orbit potential V\ (r) 
contains the long-ranged nonperturbative component, all other potentials are only relevant to the 
short range as explained by one-gluon exchange interaction. However, the observed spin-dependent 
potentials from even the latest studies [|], [7|] suffer from large numerical errors, which obscure the 
behaviors already at intermediate distances. For the phenomenological use of these potentials, it is 
clearly important to determine the form of the potentials as accurately as possible. 

For this purpose, we employ the multi-level algorithm [||] with a certain modification as ap- 
plied to the measurement of the electric-flux profile between static charges [g]. The problem is 
quite similar to this, since we need to measure the correlation function between the quark source 
and the field strength operator. This algorithm also allows us to use the Polyakov loop correlation 
function (PLCF: a pair of Polyakov loops P separated by a distance r) as the quark source instead of 
the Wilson loop. We use the field strength operator defined by F^ v = (U^y — ul v )/2i, where U^y 
is the plaquette variable. The electric and magnetic fields are then E\ = F4; and = Fjj. Noting 
((FnvFpo)) = (F^ v F pa ) pt p I '(P^ P) , where {F^ v F pa )p\ P is the two field strength correlator with the 
PLCF background, the spin-dependent potentials with the PLCF, for r = (r,0,0), are expressed as 

p 00 

V[{r) =2 dr z((B y (r,0)E z (T,z))) , (1.2) 
Jo 

poo 

V 2 \r)=2 dz z((B y (0,0)E z (T,z))) , (1.3) 
Jo 

poo 

V 3 (r) =2 dz [((B x (0,0)B x (r,z))) - ((B Y (0,0)B v (r,z)))] , (1.4) 
Jo 

poo 

V A (r) =2 dz [((B x (d,0)B x (r,z)))+2((B y (d,0)B Y (r,z)))] . (1.5) 
Jo 

It is expected from these expressions that in contrast to the use of the Wilson loop as commonly 
applied in previous works, the PLCF helps to reduce a systematic error associated with the limiting 
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Figure 1: How to construct (B y (r,0)E z (r,x)) p ip. [■ ■ ■] denotes the sublattice average. Other correlation 
functions are constructed in a similar way. 

procedure in the integral, T — » °°. At least we will have data up to T = T /2, where T is the temporal 
size of the lattice volume. 



2. Numerical procedures 

We describe the procedure how to compute the field strength correlator with the PLCF back- 
ground using the multi-level algorithm (here we restrict to the lowest level). The standard Wilson 
action is most preferable for the multi-level algorithm because its action density is locally defined. 
Thus we shall use this action in our simulation. Periodic boundary conditions are imposed in all 
directions. The essence of the multi-level algorithm is to construct the desired correlation function 
from the "sublattice average" of its components. In our case the corresponding parts are the two- 
link correlator and the field-strength-inserted two-link correlator. For a schematic understanding, 
see Fig. [j], which illustrates the computation of the correlation function, (B y E z ), for V[. 

The sublattice is defined by dividing the lattice volume into several layers along the time di- 
rection and thus a sublattice consists of a certain number of time slices. We then take averages of 
the components of the correlation function at each sublattice by updating the gauge field (with a 
mixture of HB/OR), while the space-like links on the boundary between sublattices remain intact 
during the update. We repeat the sublattice update until we obtain stable signals for the compo- 
nents. Then, we multiply these averaged components in a desired way and complete the correlation 
function. This is how the correlation function is constructed from "one" configuration. We then 
update the whole links without specifying any layers to obtain another independent gauge configu- 
ration and start the above sublattice averaging for the next configuration. 

In order to benefit from this algorithm, we need to optimize the number of time slices in 
each sublattice Afci, and the number of the internal update -/Vi up d for sublattice averaging. They 
depend on the coupling /3 and on the distances to be investigated. In principle, these parameters 
can be determined by looking at the behavior of the correlation function as a function of -/Vi up d 
for several iV ts i. An empirical observation shows that <x/V ts i = 0.3 — 0.4 fm is optimal in order to 
suppress the fluctuation of the correlation function among configurations. 
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Figure 2: Correlation functions as a function of T at r/a — 4 on the 20 4 lattice 



3. Numerical results 



We present the result obtained at jS = 6.0 (a « 0.09 fm) on the 20 4 lattice, where the ranges of 
the measured distances between static charges are R = 2 — 7. At /3 = 6.0, we found that Af ts i = 4 is 
the optimal choice. We then chose Ni up d = 7000 to be able to see the signal at least up to r/a = 7. 
The number of configuration is N con f = 76. One Monte Carlo update consists of 1 HB/5 OR. In 
Fig. |2|, we show the typical behavior of the correlation functions as a function of z at r/a = 4. 
We observe clean signals for the whole range of z. We also obtained similar clean data for other 
distances. 

Once the correlation functions are obtained, our next task is to perform the integration in 
Eqs. (1.2)— (|L5|) to obtain the potentials. Since the integration range of z is limited at most to T/2, 
we need an extrapolation to exttact the value corresponding to z —> °°. This procedure is in fact the 
potential source of the systematic error and needs careful analysis, in particular, when the statistical 
errors are significantly small as shown in Fig. f| 

Currently we applied the following analysis. Since the correlation functions were reasonably 
smooth, we firstly performed the cubic spline interpolation of the integrand and secondly evaluated 
the integral analytically in the range T€ [0, T max ], where T max = l,2,...,r/2. Then, we fitted this 
result with a function which has an asymptotic constant value at T max — > °°, like c + c'exp(— c"z max ) 
(exponential type) or c/(l + (c'/T max ) c ' ') (Hill type). The validity of the fit and the choice of the 
fitting function were monitored by looking at the the minimum of x 2 defined with the covariance 
matrix so as to take into account the correlation among different t's. The errors are evaluated from 
the distribution of the jackknife samples of the fitting parameters. 

We found that this method at least works well to extract values at T max — > °° 

forV/, V 2 'and V 3 . 

The results are shown in Figs. |3| and f| (left). For V4, however, we found that this method needs 
to be modified especially at intermediate distances, because we observed a peculiar finite z effect 
due to the symmetric behaviors of ((B X B X )) and ((B y B y }) at z = T/2 = 10. Thus, we just plot the 
integration result at Vax = 9 in Fig. ^| (right). Systematic effects of the extrapolation as well as 
finite volume effects will be investigated in future work. 
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Figure 3: The spin-orbit potentials V[ (left) and V 2 ' (right). 
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Figure 4: The spin-spin potentials V3 (left) and V4 (right). 



The qualitative behaviors of these potentials are that the spin-orbit potential V[ contains the 
long-ranged nonperturbative component (V((r) behaves as a constant at large r), while V3 and V4 
seem to be relevant at short distances. These findings are in agreement with previous works. How- 
ever, the statistical errors are significantly reduced. It is interesting to find that V 2 ' is not restricted 
to the short range, rather it has a finite tail up to intermediate distances. 



4. Summary and outlook 

We have measured the spin-dependent potentials in SU(3) lattice gauge theory with the Polyakov 
loop correlation function (PLCF) by applying the multi-level algorithm. The method presented here 
is promising to carry out further systematic investigations, such as the computation of the renor- 
malization factor of the field strength operators, Z# and Ze, as well as the scaling study, which are 
both necessary for the discussion of the fine/hyper-fine structure of the heavy quarkonium spectra. 

The preliminary studies of the renormalization factors defined a la Huntley and Michael [0], 
but using the PLCF, show the similar values as in ref. We will report these issues in our 
forthcoming publication. It is also interesting to examine the Gromes relation [pTj|], Vq = V 2 ' — V[, 
with high precision. 
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Finally we note that this method is also applicable to measuring the momentum-dependent 
potentials up to 0{\/m 2 ) JTT| , 12] from the PLCF, since the correlation functions to be measured 
are quite similar to that for the spin-dependent potentials. This will help to refine the data reported 
in ref. [^], which is in progress. 
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